Astronomy & Astrophysics manuscript no. seds ' astroph 


©ESO 2010 


September 13, 2010 





o 

Oh 
C/3 



< 

O 

Oh 
6 



(N 
> 



The RMS Survey: The Bolometric Fluxes and Luminosity 
Distributions of Young Massive Stars.* 

J. C. Mottram 1 ' 2 **, M. G. Hoare 1 , J. S. Urquhart 3 , S. L. Lumsden 1 , R. D. Oudmaijer 1 , T. P. Robitaille 4 , 

T. J. T. Moore 5 , B. Davies 61 , and J. Stead 1 

1 School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK 

2 Department of Physics and Astronomy, University of Exeter, Exeter, Devon, EX4 4QL, UK 

3 Australia Telescope National Facility, CSIRO Astronomy and Space Science, Sydney, NSW 2052, Australia 

4 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA 

5 Astrophysics Research Institute, Liverpool John Moores University, Twelve Quays House, Egerton Wharf, Birkenhead, CH41 1LD, 
UK 

6 Rochester Institute of Technology, 54 Lomb Memorial Drive, Rochester, NY 14623, USA 
Received XXX 2010 / Accepted XXX 2010 

ABSTRACT 

Context. The Red MSX Source (RMS) survey is returning a large sample of massive young stellar objects (MYSOs) and ultra-compact 
(UC) H ii regions using follow-up observations of colour-selected candidates from the MSX point source catalogue. 
Aims. To obtain the bolometric fluxes and, using kinematic distance information, the luminosities for young RMS sources with far- 
infrared fluxes. 

Methods. We use a model spectral energy distribution (SED) fitter to obtain the bolometric flux for our sources, given flux data from 
our work and the literature. The inputs to the model fitter were optimised by a series of investigations designed to reveal the effect 
varying these inputs had on the resulting bolometric flux. Kinematic distances derived from molecular line observations were then 
used to calculate the luminosity of each source. 

Results. Bolometric fluxes are obtained for 1173 young RMS sources, of which 1069 have uniquely constrained kinematic distances 
and good SED fits. A comparison of the bolometric fluxes obtained using SED fitting with trapezium rule integration and two com- 
ponent greybody fits was also undertaken, and showed that both produce considerable scatter compared to the method used here. 
Conclusions. The bolometric flux results allowed us to obtain the luminosity distributions of YSOs and UCHu regions in the RMS 
sample, which we find to be different. We also find that there are few MYSOs with L > 10 5 L G , despite finding many MYSOs with 
10 4 L G > L > 10 5 L o . 
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1 . Introduction fr om the MSX point sou rce catalogue dEgan et al.L 1 19991 120031) 

by iLumsden e t al. (2002). We define an MYSO as a radio-quiet 



. ■ Though they account for a tiny fraction of the galactic stellar mid . infrared (IR) point source which has reached a luminos . 

0> population, massive stars (M > 8 M ) dominate the galaxies . dose tQ ^ final feut where accretion is probably 

g . within which they form and evolve The radiation and outflows m { afld an ionised Hn ion hag tQ form The 

from young massive stars act as feedback processes, altering mid _ TR CQ]oms of jyrysOs are similar to those of both young 
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subsequent generations of nearby star formation. Eventually, the (tjqjj „ regions, low-mass YSOs, evolved stars, proto-planetary 
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stellar winds from these stars can lead to the destruction of the 
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servations are required to identify contaminants and gain infor- 

molecular cloud itself, ending further star formation altogether. mation about ±e candidates . The RMS survey has undertaken 

However our understanding of these effects is still in its infancy, „ j „ radio continu um observation s t0 identify radio loud UCH n 

with further study limited by the lack of large well-selected sam- ions and pNe j Urauhart et ali ftpm! [20091). ~ 1" mid-IR 

pies from which to derive the genera properties of young mas- imaging outs i de the GL IMPSE survey region faeniaminltdl 

sive stars and to study how these evolve in terms of both stellar ^| |g hurch ; ell et all Eooft to search for candidates which 

mnnn nn/j o (Tf^ 

„ . s ' , „„,, „ /r^»»oN show multiple and/o r extended sources at higher resolution 

10 th 's er^ the Red MSX Sou rce (RMS) survey dMottram et al.L [2007h and near-IR low-resolution spectroscopy 

dHoareetali |2005j | Urquhart et al| |2008b|) has carried out a se- fa order tQ id entifv evolved stars and weaki unreS ol ved UCH n 

nes of follow-up observations of a sample of candidate mas- ions ( | C larke et al.L 12003) . We have also undertaken 

sive young stellar objects (MYSOs) which was colour-selected pointed i3 CO J= j _ Q Qr J=2 _ j observations to obtain kinem atic 

* t ui m j pn i i ui i . t . ,u velocities of our sources in order to calculate kinematic distances 

Tables UJ and |2J are only available in electronic torm at the r-j -i i ^ i i 

CDS via anonymous ftp to cdsarc.u-strasbg.f r (130.79.125.5) or via flUrqunart et alj, | 2007b|, 2008a| ). 

http://cdsweb.u-strasbg.fr/cgi-bin/qcat 7J/A+A/I In a previous paper dMottram et al.L l20loh we presented new 

** E-mail:joe@astro.ex.ac.uk measurements of the fluxes of RMS sources in the far-IR (by 
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which we mean 30 /mi to 300 //m) using IRAS Galaxy Atlas 
(IGA ICao et all 1 19971) and Spitzer MIPSGAL (ICarev et all 
l2009h images. These are important if the luminosities of RMS 
sources are to be obtained, as the spectral energy distribution 
(SED) of young massive stars in both the MYSO and UCHn 
region phases peak between 60 /mi and 120 /mi. 

In this paper we will use the far-IR data from lMottram et afl 
along with other data from RMS observations and the lit- 
erature to obtain the SEDs of thes e sources ( j2] i. In {j3]we discuss 
tests of the SED model fitter of Robitaille etalJ d2007l) . which 
are used to inform the choice of input parameters for the fitter 
and evaluate what represents a good result. The results of the 
SED fitting and derived bolometric fluxes are presented in §01 
along with the calculation of filter flux ratios and flux estimates 
for those sources without far-IR fluxes derived from these ratios. 
In §|5]we present a comparison with other methods for obtaining 
the bolometric flux from SEDs, then discuss the luminosity dis- 
tributions of RMS sources in terms of their source type, before 
we reach our conclusions in §|6] 

2. Input data & the model fitter 

The SEDs of the candidate MYSOs presented in this paper are 
primarily bas ed on fluxes from v 2.3 of the MSX point source cat- 
alog ue (PSC. lEgan et all 20031) used to define th e RMS sample 
(see | Lumsden et alll2002 ) and far-IR fluxes from lMottram et ail 
( 1 201 Oh. Far-IR fluxes obtained from 18" resolution MIPSGAL 



Carey et ail l2009h 70 /mi data were used where available in 
preference to th ose from 1' x 1.7' IRAS Galaxy Atlas (IGA, 
ICaoetallfl997h data. Sources without full detections in either 
far-IR data are not included, as there is therefore no data to con- 
strain the far-IR peak of the SED for these objects (the result of 
fitting the SED without far-IR data will be discussed in 33.41 ). 
While the original coordinates for all sources were those from 
the MSX PSC, in some cases these have been updated to those 
from higher-resolution data. 

The subsections below discuss additional data included, 
where possible, to better constrain the SEDs ( {32.11 - 32.31 ), and 
provide distance information ( 32.41 ). 

2.1. Mid-1 R imaging 

In order to supplement these data, we also include fluxes 
from our ~ 1 " res olution 10.4//m TIMMI2 observations 
dMottram et all 120071) . The combined 3cr uncertainty in the off- 
set between MSX and TIMMI2 positions for all sources ob- 
served with TIMMI2 was found to be 7.6". Therefore, where 
TIMMI2 observations show no source within a radius of 7.6" of 
the source coordinates, a non-detection upper limit flux is used 
based on the uncertainty in the 10.4/mi TIMMI2 flux. 

From TIMMI2 and other high-resolution mid-IR imaging, 
it is clear that some of the flux detected within the MSX beam 
(18") in star formation regions is due to diffuse emission not 
directly relating to the source of interest. Therefore, the bolo- 
metric flux of larger-beam observations should be apportioned 
based on the higher-resolution data. In these cases, the apportion 
ratio was therefore calculated as the ratio of the 10.4 /mi flux 
to the MSX 8 //m flux. Where the TIMMI2 flux is larger than 
the MSX 8 fim flux, the ratio is set to 1 . The individual fluxes 
from large-beam observations (i.e. MSX, MIPSGAL, IGA, sub- 
mm and mm) were then multiplied by the apportion ratio before 
input into the model fitter and the uncertainties of these fluxes 
calculated as the addition in quadrature of the uncertainty in the 
apportion ratio and the measurement uncertainties. 



In using a simple ratio of the 10.4 /mi flux to the MSX 8 /mi 
flux we assume that all the sources detected in the TIMMI2 ob- 
servations have the same SED shape and that the 9.7 /mi sili- 
cate feature affects both fluxes equally. While this is an obvious 
simplification, it is certainly an improvement over not perform- 
ing this correction and the absence of higher resolution data at 
longer wavelengths for all our sources makes a more detailed 
determination impossible at this time. This method was required 
for 178 out of a total of 1 183 sources. 

Similarly, some RMS sources have detections at 8 //m in the 
GLIMPSE PSC, in which case this flux is used to calculate the 
apportion ratio as for TIMMI2 sources. 182 out of a total of 
1 183 sources were apportioned in this way. In addition, in some 
cases it is clear from higher-resolution mid-IR imaging that more 
than one YSO or H n region contributes significantly to the flux 
within the MSX beam. For these sources, the individual compo- 
nents are given letter designations and apportioning carried out 
where possible. Where we do not have measured fluxes, such 
as cases which can be seen in the GLIMPSE 8 /mi images but 
do not have 8 //m fluxes in the PSC, it is still clear that each in- 
dividual source does not contain the total MSX flux. Therefore 
the remaining fraction of the total flux of the larger-beam obser- 
vations is split equally between those sources which cannot be 
apportioned using TIMMI2 or GLIMPSE 8 /mi fluxes. This is 
performed for 136 out of a total of 1 183 sources. 

2.2. 2MASS 

In order to help constrain the near-infrared side of the SED, 
J, H and K s band fluxes from the nearest 2MASS PSC 
dSkrutskie et all 120061) entry within 7.6" (see above) are used. 
Though this radius may seem large, it is still less than half the 
full-width half maximum (FWHM) of MSX (18") and is moti- 
vated by the uncertainty in the source coordinates. In the case 
of sources w ith TIMMI2 observat ions, the 2MASS counterparts 
identified bv lMottram et al.l d2007l) are used instead. Where there 
is no 2MASS detection within 7.6" of the MSX PSC coordi- 
nates, or where inspection of the 2MASS images during object 
classification indicated that the nearest source is unlikely to be 
the near-IR counterpart to the MSX target, the noise in each band 
is calculated in the following way in order to estimate an up- 
per limit flux. All 2MASS PSC entries which are detected in 
other bands but not in the band in question within 1' of the target 
coordinates are identified. An average of the upper limit fluxes 
in the band in question of all such sources is then calculated. 
While some 2MASS counterparts may be mis-identified by this 
method, the SEDs of young massive sources are dominated by 
emission in the far-IR, so this is unlikely to have an overriding 
influence on the calculated bolometric flux. 



2.3. Sub-millimetre and Millimetre Continuum 

SCUBA 450 //m a n d 850 //m data from the legacy archive 
dDi Francesco et all 120081) are included, where available, for 
those sources with an offset between the SCUBA observations 
and the source coordinates < 14", the beam full-width at half 
maximum (FWHM) of SCUBA at 850/mi. The FWHM is used 
rather than the pointing error (usually of order + 5") as many 
of these SCUBA sources associated with RMS young massive 
stars are larger than the beam size. In cases where the offset 
was larger than 14" but smaller than the source size at 850/mi, 
the SCUBA fluxes were used as upper limits to the true flux 
at these wavelengths. The percentage errors on the SCUBA 
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450 Aim and 850 //m fluxe s are 50 
(iDi Francesco et al.Ll2008h . 



and 20 % respectively 3. Fitting the spectral energy distributions of young 

massive stars 



Millimetre continuum observations at 1 .2 mm with the SEST 
imaging Bolometer array (SI MBA) on the Swedish ESQ sub- 
millimetre telescope (SEST) by IFaundez et all d2004l) . lHill et al l 
d2005l) and lBeltran etail fcOOd) . along with 1 .2 mm observations 
with the Max-Planck Millimetre Bolometer (MAMBO) on the 
Institut de Radio Astronomie M illimetrique (IRAM) 30 m tele- 
scope by IBeuther et al.l d2002l) were also included. Originally, 
all sources with offsets from the source coordinates less than or 
equal to the beam size (24" for SIMBA and 11" for MAMBO) 
were included as detections and those with offsets larger than 
the beam but smaller than the source size as upper limits based 
on the source flux. However, in some cases the source sizes as 
derived by Gaussian fits are large compared to the beam size 
and/or the source is irregular and extended at 1.2 mm, resulting 
in a much higher flux at this wavelength than is consistent with 
the rest of the data. All 1 .2 mm images were therefore examined 
and compared to MSX images in order to determine which ob- 
servations should be regarded as upper limits to the true source 
flux. 

The uncertainties on detected sources are ta ken as the rms 
noise from the vario us data sets, i. e. 40 mJy IFaundez et alj 
2004. 40 mJy fo r IBeltran et al.l (l200j 



(200 4J), 40 mJy fo r IBeltran et al.l (120061) and 1 5 mJ y for 
IBeuther et all (120021) . The observations of lHill et all (120051) were 
observed over multiple runs, so in order to be conservative the 
largest rms noise value from these runs (i.e. 150 mJy) was used 
for all data. For non-detections, or sources with offsets from the 
MSX position larger than the measured size of the detection at 
1 .2 mm, 3 cr of the appropriate minimum background flux levels 
are used as an upper limit to the 1.2 mm source flux. 



2.4. Kinematic Distances 

Kinematic distances have been obtained for the majority of 
RMS sources, calcul ated from kinematic velocities using the 
iBrand&Blitzl (119931) rotation curve, assuming the galactocen- 
tric distance and velocity of the solar system to be 8.5 kpc 
and 220 kms 1 respectively dUrquhart et all l2007bl l2008ah . All 
sources with a galactocentric radius < 8.5 kpc have a kine- 
matic distance ambiguity caused by the fact that two points 
along the line of sight, one near and one far, have the same 
line of sight velocity in the galactic rotation curve, though 
at the tangent point these are the same and so not an is- 
sue. We have s olved this distance a mbiguity in the north- 
ern hemisphere dUrquhart et all 1201 Ol) by cro ss-matching our 
source s with a sample of clouds identified by iRathborne et al.l 
d2009l) in the 13 CO Boston University Five College Radio 
Astronomy Observatory (BU-F CRAO) Galactic Ri ng Survey 
(GRS, 18° < 1 < 55.7°, |b| < 1 °: Uackson etail l2006h for which 
the am biguity has already been solved by Roman-Duval et all 
d2009l) . In the southern hemisphere we are using data from the 
South ern Galactic Plane Survey (SGPS. lMcClure-Griffiths et all 
120051) and the H i self-absorption method (ILiszt et all 119811) to 
solve the ambiguity to our sources (Urquhart et al, 2010b, in 
prep.). 

For sources where there are unsolved near-far distance am- 
biguities, fits are calculated for all distances. These 33 sources 
are included in the Table of results (table [TJ presented in §|4] as 
the average of these fits, but are not included in the analysis. 



This section discusses how the bolometric flux is obtained from 
model fits to spectral energy distributions ( 33. Il l, then describes 
a test designed to investigate how use of an incorrect distance 
in the fitting process affects the accuracy of the resulting bolo- 
metric flux ( 33.2l i. The choice of the Ay range that the fitter is 
allowed is discussed ( 33.31 l. followed by a study of the relative 
importance of the various data sets that make up the SED data 
( 33.4I >. The investigations in 33.21 and 33.41 were run using the 
data for G034.7569+00.0247 as the data set for this probable 
MYSO is one of the most complete. 

3.1. The model SED fitter 



The radiative transfer models of iRobitaille et al.l d2006l) cover a 
range of different stellar, disk and envelope prop erties at a range 
of ang les of inclination. The model fitter of IRobitaille et al.l 
(120071) then takes these models, convolves them with a user- 
specified set of filters, allows a range in distance and Ay and then 
minimises the x 2 between the model filter fluxes and the input 
data. We use this model fitter in preference to other methods as it 
takes into account dust features, such as the 9.7/iin silicate fea- 
ture, and consider a more realistic geometry for the circumstellar 
material than a simple sphere. It is important to mention that a 
good fit to the data does not prove that the geometry or condi- 
tions present in the model are a good description of the observed 
source, rather that the model is consistent with the observational 
data. In particular, it may be the case that the data are insufficient 
to fully constrain the 14 parameters used to construct the models, 
so it is not appropriate to infer knowledge of variables t hat th e 
SEP data alone do not specifically constrain dRobitailld, 120081) . 
Ide Wit et all d2010l) found that the best fitting models returned 
by the fitter for MYSOs usually have quite massive dust disks 
which are ruled out by mid-IR interferometric observations. In 
addition, due to the non-random nature in which the parame- 
ters were sampled care must be taken to establish that observed 
tren ds in the data are not a result of trends in the parameter space 
(see lRobitaiiieLl2008L for a full discussion). 

The bolometric flux of the fit is, howev er, well constrain ed by 
the fits, as this relates directly to the SED dRobitailld. l2008h . The 
following analysis will therefore focus on this property, while 
comparisons between the bolometric flux obtained by SED fits 
with other methods will be addressed in 35.11 

It is important that values derived from the results of the 
model fitter include several models, in order to obtain a reason- 
able estimate of whether these results are constrained and unique 
with the data used. Therefore the number of fits within a Ax 2 =2 
of the best fit was calcul ated (we use the same definition of x 1 
as Robitai lTe et alll20 07l). If there were less than 10 fits included 
within x\ st + AAf 2 tnen AA' 2 was increased by integer steps until 
this was no longer the case. This method was used rather than 
simply finding the best 10 fits so that the spread in fits was not 
artificially limited if many models are good fits to the data with 
a range of values for distance, bolometric flux and Ay. A mini- 
mum of 10 fits were required in order to ensure that one or two 
fits did not strongly bias the results. 

The model bolometric flux for the source was obtained by 
using the model luminosity and model distance derived for each 
individual fit to calculate the bolometric flux for each fit (/,). 
The luminosity returned by the model fitter is the combination 
of that for the central source, accretion onto the central source 
and accretion in the disk, which is then corrected for foreground 
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A (/Lim) 



A (/xm) 



Fig. 1. Example SED fits with allowed distances of 2.0 ±1.0 kpc (left) and 5.0 ±1.0 kpc (right) respectively. The circular points 
indicate flux detections while the downward pointing triangles indicate upper limit fluxes. The solid line shows the best fit to the 
data, while the grey lines indicate other fits within a given Ay 2 of the best fit. 
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Fig. 2. Bolometric flux (left) and Ay (right) vs. the mean model distance obtained from SED fits with specified model distances 
d ± 1.0 kpc with d varying between 2 kpc and 13 kpc for the data of G034. 7569+00. 0247. The black horizontal error bars indicate 
the uncertainty on the mean model distance while the dashed horizontal error bars indicate the range of allowed model distances. 



extinction consistent with the SED, and is therefore independent 
of viewing angle unlike the SED itself. Next the weighted mean 
and weighted standard deviation were calculated of all fits within 
x\ st + Ay 2 to give the flux and flux uncertainty using: 
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where N is the number of non-zero weights, and the weights are 
calculated from the x 2 as: 
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The x 1 values are used to weight the fits so that the best model 
fitted to the data has more impact on the results than poorer fits, 
but that the results are not heavily skewed away from the mean. 
This method is also used to calculate the mean model distance, 
the mean Ay, and mean filter flux ratios given the results for 
individual model fits to the data. 



3.2. Model Distance 

A distance range is required by the model fitter within which 
to fit the observed data. Therefore the test described below was 
undertaken in order to explore the distance dependence of the 
bolometric flux obtained from the model fitter, and thus what 
effect an incorrect distance has on the model derived flux. 

Fits were performed to the data for G034.7569+00.0247 for 
model distances d + 6d kpc, with d varying between 2 kpc 
and 13 kpc in 1 kpc integer steps. Example SED fits for 
2.0 ±1.0 kpc and 5.0 ±1.0 kpc ( the kinematic distance for 
G034.7569+00.0247is 5.0 ± 0.7 kpc. lUrquhartetaill2008al) are 
shown in Figure[T](left and right respectively). For all runs of the 
model fitter, 5d is set to 1.0 kpc as this is the order of the typical 
error in the kinematic distance for RMS sources. 

The model fitter divides the user-specified distance range 
logarithmically in order to obtain model distances after which 
the model SED fluxes are convolved with common filter bands 
and interpolated to the apertures used to obtain the data. The 
logarithmic binning was set to 0.01 for this test, while 0.025 was 
used for the other tests and the results presented in §|4] The re- 
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100 



Fig. 3. Left: A histogram of Ay derived by the model fitter for all 
RMS sources with 6F W / F w < 0.5 andx 2 best < 30 for Ay 2 < 10. 



suits were checked for a few sources at the lower model dis- 
tance binning step and the differences were found to be min- 
imal. At small model distances the apertures may not encom- 
pass the whole model but only the inner region, the calcula- 
tion of which can require large amounts of computer memory. 
Therefore, where the distance to the source was less than 2.4 kpc, 
the distance given to the fitter for the main results was reset to 
this minimum limit. 

The results of the fits where d is varied for the data of 
G034.7569+00.0247 are shown in Figure |2] in terms of bolo- 
metric flux (left) and Ay (right) vs. the mean model distance as 
calculated using the weighted mean and weighted standard de- 
viation of the model distances for the individual fits within A^ 2 
of^? ,. The Ay is relatively constant with distance for these fits 
(see Figure [2]), suggesting it is constrained by the data. The in- 
crease in the bolometric flux and error at small distances is due 
to the source being so close that the Ay is insufficient to stop 
emission in the visual being observed (see Figure[TJ. The model 
fitter allows this as the near-IR fluxes are upper limits, and so are 
treated as weaker constraints on the SED than the detections at 
longer wavelengths. Overall, though the bolometric flux appears 
to decrease with increasing distance, the change is of similar or- 
der or smaller than the uncertainties in the fluxes, and so the 
results can be considered distance independent. 
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3.3. Visual Extinction 

The visual extinction (Ay) between Earth and the edge of the cir- 
cumstellar material considered by the model is an important free 
parameter in fitting model SEDs to observed data. Restricting 
this value too much could result in poor fits, particularly be- 
tween the near/mid-IR and far-IR/sub-mm parts of the SED, 
while allowing too much freedom could lead to unphysical so- 
lutions. The extinction law model assumed was calcula ted by 
iRobitaille et all d2007l) using the method of iKim et alj (fl994l) 
for a typical galactic IS M, taking into account the results of 
llndebetouw et al] (120051) with regards to mid-IR extinction. 

We limit the Ay to be within the ra nge to 50, in lin e 
with typical values for UCH n regions (e.g. lHanson et al.ll2002l) . 
though as shown in Figure [3] the majority of sources have fitted 
values below this. 



io- ,4 r \ 

0.1 1 10 100 1000 10* 

A (/im) 

Fig. 4. SED fits with different data sets included. Top: All data 
included. Middle: All data except MIPSGAL. Bottom: All data 
except sub-mm and mm. Bottom right: All data except 2MASS. 
The lines and points have the same meaning as in Figure Q] The 
bolometric flux derived for each data set is indicated on the plot. 



3.4. Data Included 

SED model fits for G034.7569+00.0247 including only subsets 
of the available data were undertaken in order to explore the ef- 
fect of missing data on the bolometric flux results, examples of 
which, along with the fit using all data are shown in Figure [4] 
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Comparing the top plot, which includes all data, with the mid- 
dle plot where the MIPSGAL data has been omitted, the effect 
of using the model fitter for sources without far-IR data is re- 
vealed. The lack of constraint at the peak of the SED leads to the 
fractional error in the bolometric flux being larger than 1, i.e. it 
is unconstrained, despite having an otherwise reasonably com- 
plete data set. As a result, as mentioned in §f2]we do not fit SEDs 
for sources where we do not have far-IR flux data, though we 
will see in 34. II that we can use the F21 to Fb i ratio as derived 
from fitted sources to estimate Fb i for these sources. The lower 
plot of Figure|4]illustrates that the lack of sub-mm and mm data 
for many of our sources does not change the results within the 
calculated errors. 

4. Results 



The SED model fitter of iRobitaille et alj d2007l) was run for all 
1183 sources with far-IR data as discussed in fusing kine- 
matic dist ances with 6d - 1 kpc obtained from molecular emis- 
sion lines dUrquhart et al.L l2007bl l2008al 1201 Ol Urquhart et ai, 
2010b, in prep.) toward each source. There are 33 sources re- 
maining with multiple solutions for the kinematic distance to 
some sources due to unresolved near/far distance ambiguity, for 
which the SED fitter was run for each possible distance and the 
flux and error are given as the mean and <x/ -\Jn of the results from 
these runs. The maximum Ay allowed was 50 (see 33.31 ). 

A sample of the results from these fits are shown in Table [1] 
while the full version of this table is available online at the 
CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.125.5) 
or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ The RMS 
sample sources have been given one of a number of clas- 
sifications, shown in column 4 of Tab le [T] which are dis- 
cussed in §2.1 of iMottram et al.l d2010t) . 'Young / Old Star' 
is abbreviated to 'Y/O Str' and an additional classification 
for more evolved Compact or diffuse Hn regions, where the 
ionised region fills much of the MSX beam, is abbreviated 
to 'CHn'. For details on individual sources see the RMS 
database (http://www.ast.leeds.ac.uk/RMS). As distance solu- 
tions are subject to change as new information becomes avail- 
able, these are not included in the table. However, up to date dis- 
tance information is available on the database and once distances 
are available for all sources, these will be the subject of a future 
publication. In addition, as discussed in 33.2l the flux results are 
distance independent. The data included in the SED of a source 
is given by a series of flags in column 6 of Table [1] These are 
from left to right 2MASS (Y or N), TIMMI2 (Y or N), whether 
the far-IR was from the IGA (I) or MIPS GAL (M), submillime - 
tr e (Y or N) and millimetr e data (N 1 =lBeufher et alJ d2002l). 
2=lFaundez et all d2004h . 3= lHill et all d2005h . 4= lBeltran et alJ 
d2006j)). Table Q] contains results for 1173 sources, while for a 
further 10 sources the fitter was run but returned a flux smaller 
than the error in the flux (i.e. unconstrained). 

We do not discuss in this paper the results o f the 14 free pa- 
rameters in the models of lRobitaille et al.l d2007h as these are not 
uniquely constrained by the data, and the models do not cover 
the parameter space uniformly. Indeed, the distributio ns of most 
of the parameters follow trends in the models (see iMottraml 
I2008L for more details). 

A total of 1069 sources have an unique assigned kinematic 
distance and have fits with 5F W / F w < 0.5 and x\ t < 30 for 
AX 2 < 10, of which 488 are identified as YSOs, 20 as H n/YSOs, 
45 as Young / Old Star, 76 as CHn and 440 as Hn regions. 
Histograms of the number of sources as a function of their 
luminosity for YSOs and H 11 regions are shown in Figure [5] 




1 000 1 o 4 

Luminosity (L©) 

Fig. 5. Plots of the number of sources vs. luminosity for RMS 
candidates with classifications of YSO (488 sources, black) and 
Hn (440 sources, red). The blue arrows in dicate the lum i nosity 
of a B0 type star from the calculations of iMartins et al.l (I2005I) 
corrected for the l umino sity/temperature/mass relationship of 
iMevnet & Maederl d2000l) . 



These sources have typical uncertainties in the bolometric flux 
of ~ 17%, have a median distance of 5.1 kpc, a typical distance 
uncertainty of ~ 19 % (assuming 6d = 1.0 kpc), and have typ- 
ical uncertainties in the luminosity of 34%. The sources with 
MIPSGAL far-IR fluxes do not show any significant difference 
in terms of distribution to those with IGA data, implying that 
the improvement in resolution provided by the MIPSGAL data 
does not result in a lack of luminous sources. Figure [6] shows 
the plots for YSOs and H 11 regions from Figure [5] for both the 
assigned distances (black) and with all near/far determinations 
for all assigned kinematic distances inverted. This is undertaken 
to examine the effect of incorrect near/far distance ambiguity as- 
signments on the luminosity distributions. As can be see, this 
does not change the overall shape of the luminosity distributions 
appreciably. Therefore the remaining individual uncertainties in 
the distances to the RMS sample do not significantly affect the 
global results on a statistical level. These results will be dis- 
cussed in detail in $5] 



4.1. Sources Without Far-IR Data 



Though far-IR fluxes were obtained by IMottram et al.l d2010h 
for the majority of young RMS sources, this was not possible 
for some sources due to confusion within the available data, 
which in turn prevented SED fitting of these sources (see 33.4l i. 
However, the properties of those sources where fits were ob- 
tained can be used to calculate the mean ratio between the bolo- 
metric flux and the MSX 21 jjm flux for all sources with good 
SED fits. This mean ratio can then be used to obtain estimates 
of the bolometric flux of RMS sources which do not have far-IR 
data, since all have measured MSX 21 /mi fluxes. 

In addition to the total flux, filter-band fluxes were output for 
the SED model fits to the data. The MSX 21 /mi band filter fluxes 
in WrrT 2 can be calculated by multiplying the SED MSX 21 /mi 
flux in Janskys by t he bandwidth of the filter (4.041 x 10~ 14 Hz, 
ICohen et all |200 0j). The weighted mean ratio of the total bolo- 
metric flux to the MSX 21 /mi band filter flux can therefore be 
obtained for each source, in a similar manner to that used for 
the bolometric flux (see §0, with the uncertainty given by the 
weighted standard deviation. The results of these calculations 



6 



J. C. Mottram et al.: The Bolometric Fluxes of Young Massive Stars 



Table 1. Example results of SED fitting for all sources with far-IR data. The columns are as follows - 1: MSX PSC name. 2: 
Right ascension. 3:Declination. 4: RMS classification. 5: Number of detected flux data points, not including upper limits, included 
in the fit. 6: Flags indicating which data were included in the SED. These are from left to right 2MASS (Y o r N ), TIMMI2 (Y 
or N) far -IR IGA (D or M IP SGAL (M), submilli metre (Y or N) and millimetre (N, 1 = lBeuther etail d2002l) . 2= lFaundez eTall 
d2004l) . 3= lHill et al.1 J2005I) . 4= lBelfran et al] d2006l) ). 7: x\ est for the fit. 8: A* 2 at which at least 10 fits are found. 9: The number 
of fits. 10: The bolometric flux with uncertainty, calculated using equations [TJ and [2] 11: Whether the kinematic distance has been 
solved (s) or the near/far ambiguity remains unsolved (nf). Sources with galactocentric radii > 8.5 kpc or at the tangent point are also 
indicated as solved (s). 12: The apportioning correction applied to large-beam fluxes, as discussed in 32.11 with a value of 1 meaning 
that no apportioning was performed. 13: Type of apportioning used with 'N' for none, 'T' for TIMMI2 data, 'G' for GLIMPSE 
data and 'S' for a simple split of the remaining fraction between all multiple source components without TIMMI2 or GLIMPSE 
detections. A full version of this Table is available online at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79. 125.5) or via 
http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ 
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Fig. 6. Plots of the number of sources vs. luminosity for YSOs (left) and Hn regions (right) with the assigned distances (colours 
as in figure[5]l and with all near/far distance determinations swapped (blue and green respectively). The blue arrows have the same 
meaning as in Figure [5] 



do not appear to be dependent on source flux or type (see left- 
hand plot of Figure[7]), though the scatter in sources is relatively 
large and so may mask such relationships. The ratio itself has 
a roughly log-normal distribution, as shown in the right-hand 
plot of Figure 13 Sources with MIPSGAL data generally have a 
slightly smaller spread than the general distribution, so the prop- 
erties of the Gaussian fits to the mean ratio for these 613 sources 
are used. 

This mean ratio was therefore used to obtain estimates 
of the bolometric flux of RMS sources which do not have 
far-IR data, since all have measured MSX 21/vm fluxes. The 



bolometric fluxes obtained for sources with TIMMI2 (69) or 
GLIMPSE PSC (20) fluxes, as well as split sources (19), 
were apportioned using the method discussed in 32.11 An 
example of these results is presented in Table [2] while the 
full version of this table is available online at the CDS via 
anonymous ftp to cdsarc.u-strasbg.fr (130.79.125.5) or via 
http://cdsweb.u-strasbg.fr/cgi-bin/qcat?I/A+A/, 

Estimate total fluxes were obtained using this method for an 
additional 280 young RMS sources. The uncertainty in the total 
flux for each source is a combination in quadrature of the uncer- 
tainty in FboI / Fmsx 21 and the uncertainty in the MSX 21 //m 
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Bolometric Flux (Wm 2 ) F Bol / F usx „ Mm 

Fig. 7. Left: The distribution of F% i / F msx 21 as a function of Fr \ for all YSOs and H11 regions (red and blue respectively in 
the online version of this Figure). Error bars are shown in black in order to give an idea of the total variation in values. Right: The 
histogram of Fb q \ / Fmsx 21 for all sources with MIPS data. A Gaussian is fitted to the data is shown by the dashed line. 
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Table 3. Filter flux ratios derived from the SED fits to sources 
with MIPSGAL detections. 



flux. The fact that the Gaussian fit to Fr \ / Fmsx 21 is performed 
in log space results in asymmetric uncertainties, which are prop- 
agated through to the total flux. 

The distributions of luminosities for sources without far-IR 
data identified as YSO or Hn region with uniquely assigned 
distances is shown in Figure [8] with the distributions of SED 
derived luminosities from Figure [5] also shown for comparison. 
Due to the lower number of sources (110 YSOs and 103 H 11 re- 
gions), the bin size for the distributions for luminosities derived 
using F B oi / Fmsx 21 is twice that used for the SED luminosities. 
The spread of sources is larger than for the SED luminosities, 
probably caused by the larger uncertainties in the bolometric 
fluxes derived using the MSX 21 yum flux and the fact that these 
sources are in regions which are confused at far-IR wavelengths. 
However the general shape of the distributions is similar, and the 
peaks in the distributions are near those for the SED luminosi- 
ties. This method of estimating the far-IR flux therefore seems 
to provide a reasonable estimate when no far-IR information is 
available. 



4.2. Filter Effects 

As with the ratio of F^oi I Fmsx 21 discussed above, various 
other ratios both of filter fluxes to the bolometric flux and fil- 
ter fluxes to other filter fluxes can be calculated from t he mo del 
SED fits. However the SED models of lRobitaille et al.l (120061) do 
not include PAH emission features, so the flux ratios for bands 
where these are important may be different, e.g. the 8.0/im PAH 
band may affect ratios involving the MSX 8 fim band. The results 
of these fits are presented in Table [3] 



5. Discussion 

5. 1 . Comparison with Other Methods 

In order to explore the relative accuracy of other methods com- 
monly used to obtain bolometric fluxes compared to the SED fit- 
ter results, we performed calculations using both simple trapez- 
ium rule integration and a combinati on of two greybody fits (e.g. 
iMinier et all l2005t iHiU et al.Ll2009l) with model rather than ob- 
served fluxes. The SED model 'observed' flux for each fit was 
first corrected for the extinction derived by that fit. Next the 
weighted mean and weighted standard deviation were calculated 
in order to obtain a corrected model flux with error for each filter 
where data were input to the fitter for that source. The greybody 
fit was obtained using two temperature components given by: 

F v = B y (TiXl - e-^)e\ + B v {T 2 ){\ - e^) 0\ 

where B V (T) = -n^—, and t v = T ref (v / v ref f ^ 

where 8 is the source size, T is the dust temperature, r v is the 
optical depth at frequency v and B is the dust emissivity index 
such that t can be evaluated at any frequency relative to a ref- 
erence frequency v re /. A total of five free parameters were used 
for the fits: the temperatures corresponding to the peaks of the 
two greybody components at lower and higher temperature (Ti 
and T2), the reference optical depth (j re f) at 850 fim, the source 
size in arcseconds of the high-temperature component (62) and 
B\. B2 was kept constant at 1, while 6\, the source size of the 
low-temperature component, was set to the beam size of 1.2 mm 
SEST data (i.e. 24"). B\ was restricted to only be between 1 and 
2. Only sources with at least five data points were fitted using this 
method, due to the number of free parameters, but this resulted 
in the exclusion of only a few sources which tend to have only 
upper limit 2MASS fluxes, an upper limit in the MSX 12/mi 
band and MIPSGAL data. Though these sources have less than 
5 detected data points, the SED fitter is also constrained by the 
upper limits. 

Figure [9] shows a comparison between bolometric fluxes de- 
rived by the model SED fitter with those obtained by simple 
trapezium rule integration (left) and two component greybody 
fits (right). There is quite a large scatter, with a significant frac- 
tion more than 50% away from the SED fitter derived flux. Both 
alternative methods are more likely to over-estimate the flux, 
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Table 2. Example bolometric fluxes for young RMS sources without far-IR data. A full version is available online at the CDS via 
anonymous ftp to cdsarc.u-strasbg.fr (130.79. 125.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/, 
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Fig. 8. Plots of the number of sources vs. luminosity calculated from MSX 21 /jm fluxes for RMS candidates without far-IR data 
with classifications of YSO (110 sources, left) and Hn region (103 sources, right). The blue arrows have the same meaning as in 
Figure|5] The dot-dashed histograms show the normalised luminosity distributions from Figure[5] 



with mean values of F grev / ¥ finer = 1-3 5 and ¥ trap /Ffn, er = 1.16. 
If the Emerson flux (lEmersonl 1 1988b is calculated with IRAS 
fluxes from the models instead of a full trapezium integration, 
the mean value of F Emerson / ¥f itter = 1.09 though this does not 
include the effect of contamination which is often a major factor 
for IRAS fluxes of sources in the galactic plane. If IRAS PSC 
fluxes are used, ¥ Emerson I F filler = 1-82 with considerably more 
scatter. Many of the greybody fits do not fit both the near and 
mid-IR well, as shown by the example in Figure [10] The group 
of sources with considerably larger greybody fluxes in FigurefTUl 
are all sources with MIPS far-IR data but no submm or mm 
points where the fit has the far-IR peak at A > 70yum, and so 
overestimates the flux from cold dust. Neither method will deal 
with the 9.7 fim silicate feature accurately. FigureQT]shows a his- 
togram of the fractional errors of fluxes derived by the SED fit- 
ter, greybody fits and trapezium integration. These indicate that 
errors in fluxes derived using trapezium integration are system- 
atically larger than those from obtained using the other methods, 
and that the errors using greybody fits are on average larger than 
those using the SED fitter. The median fractional error is 0.17 
for the SED fitter fluxes, 0.24 for the two-component greybody 



fluxes and 0.59 for trapezium integration. Overall, therefore, the 
SED fitter produces more consistent and accurate results than 
either greybody fitting of simple integration. 



5.2. YSOs 

The upper left plot of Figure [5] shows that the number of YSO 
sources peaks around 5 x 10 3 L Q , and decreases sharply above 
~ 2 x 10 4 L . The gradual drop in the number of YSOs at 
luminosities below the peak is caused by incompleteness in 
the RMS sample due to the 21 yum flux-limit of the MSX sur- 
vey. Though there are some sources which have L > 10 L Q 
(which corresponds to an ~ 30 M Q 7 type star using the tem- 
perature/logQ scale of iMartins et al.1 d2005l) corrected for the 
lumino sity/temperature/mass relationship of iMevnet & Maederl 
j2000h ). inspection reveals that many of these sources should be 
apportioned based on GLIMPSE images, and thus have lower 
luminosities, but lack PSC detections at 8 ^m. As already dis- 
cussed in §|4] the shape of the overall distribution is robust 
against these issues. 
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Fig. 9. Comparison of bolometric fluxes derived by the model SED fitter with those obtained by simple trapezium rule integration 
(left) and two component greybody fits (right). The green dashed line indicates y = x. The blue dashed lines indicate equality ± 50%. 
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Fig. 11. Histogram of the fractional errors (i.e. 6Fb \ / Fb i) for 
fluxes obtained using the SED fitter (black), trapezium integra- 
tion (red) and greybody fitting (blue). 

5.3. Hn Regions 

The left-hand plot of Figure [5] shows that the number of H n re- 
gions in the RMS sample peaks around 3 x 10 4 L Q . There are few 



sources which are identified as H n regions but have a luminos- 
ity derived from the SED lower than the lim it for formation of an 
H n region ( L ~ 550 L , s pectral type ~B3. lMartins et alll2005b 
iMevnet & Maederl 120001) . Inspection of these sources suggest 
that they are correctly identified as H n regions, but suffer from 
individual distance determination issues. 

While the models used by the SED fitter do not include ad- 
ditional sources of heating of th e dust inside the io nised zone, 
in particular Lyman alpha (e.g. iHoare et aU 119911) . the SEDs 
of UCH ii regions are similar to those of YSOs except at radio 
wavelengths, which the models do not probe. As the primary 
concern of this work is the bolometric flux, the fits are accept- 
able for this purpose. 

The most luminous UCHn region observed by 
IWood & Churchwelll (Il989t) has a luminosity of ~3 x 10 6 L , 
while no R MS source has a luminosities above 10 6 L . However, 
some of the lWood & Churchwelll d 19891) sources are now known 
to contain multiple objects and in some cases the RMS kine- 
matic distance measurements are si gnificantly less than the ones 
used by those authors. For example lWood & Churchwelll (1 19891) 
use d = 9.7 kpc for G045 .4543 +00.0600 but the RMS kinematic 
distance to this source is 7.3 kpc, which would change their 
1.44 x 10 6 L to 8.2 x 10 5 L . The SED fit for this source 
gives L = (3.7 + 1.0) x 10 5 L but MIPSGAL rather than IRAS 
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Fig. 12. SED (left) and MIPSG AL 70 urn image (righO for t he Hn region G045.4543+00.0600, for which we derive a luminosity 
of (3.7 + 1.0) x 10 5 L while I Wood & Churchwelj (i 19891) derive a luminosity, after correction to our distance of 7.3 kpc, of 
8.2 x 10 5 L using IRAS far-IRdata. The RMS source is merged with the source to the east (IRAS 19120+1103) in IRAS data and 
has a PSC 60 fim flux of 5340 Jy compared to a MIPS 10 fim flux of 2402 Jy, accounting for much of the remaining difference in 
luminosity measurements. Saturated pixels are indicated in yellow in the online version of this Figure. For the SED, the lines and 
points have the same meaning as in Figure Q] 



far-IR fluxes are used and the MIPSGAL 70 fim image shows 
multiple bright regions nearby (see Figure fT2l . The SED does 
not fit the 850 fim point well, but this source exhibits significant 
surrounding low-level emission in the MIPSGAL image and 
strong radio emission. Both of these factors could increase 
the observed SCUBA flux above that consistent with the SED 
models based on the lower wavelength emission, which is well 
fit. 



5.4. Comparison of YSOs and Hn Regions 

In comparing the luminosity distributions of YSOs and H n re- 
gions (Figure [5}, the difference in the peak of these distribu- 
tions is apparent. Indeed, a Kolmogorov-Smirnov (KS) test of 
the cumulative distribution functions, indicates that the proba- 
bility of these distributions being the same is of the order 10~ 26 . 
It is unlikely that this is due to an under-detection of B spectral 
type H n regions in the samp le as the RMS radio observations 
(lUrauhart et all 120071 12009T) were deep enough to detect Hn 
regions down to spectral type Bl to distances of 20 kpc assum- 
ing an electron temperature of 10 4 K. In addition, many of these 
MYSOs are visible at near-IR wavelengths, so are unlikely to be 
heavily embedded hyper-compact (HC) H n regions. 

Main-sequence stars with a spectral type of B3 (L > 550L o ) 
or earlier are capable of producing H n regions, and evolutionary 
tracks indicate that core hydrogen burning should have started in 
YSOs of spectral type B3 or earlier despite t he fac t that accretion 
is still ongoing (e.g. lYorke & Sonnhalterl 120021) . It has there- 
fore been unclear until recently why MYSOs, where the cen- 
tral source is massive enough to form an H n region but no H n 
region is observed, are detected at all, rather than all young mas- 
sive stars transitioning directly from i ntermediate mas s YSO s 
to UCHn regions as they gain mass dHoare & Francol 120071) . 
Some mechanism must either qu e nch o r mask the e xistence 
of an Hn region (e.g. IWalmslevL 119951: iKetol 120031) . or in- 
crease the radius of the star such that the stellar T e ff, and 



thus UV flux, are decreased (e.g. Yorke & Bodenheimerl |2008; 



IHosokawa & OmukaiL l2509t IHosokawa et all 120 1 oTT 

Ascertaining further how selection effects influence the dif- 
ference in the distribution of RMS YSOs and UCH n regions will 
require the luminosity functions of these sources, which will be 
the subject of a future paper (Mottram et al., 2010, in prep.). 

6. Conclusions 

After preforming tests of the influence of the i nput parameters to 
the model SED fitter of lRobitaille et alJd2007l) . we have obtained 
the luminosities of 1173 young RMS sources, of which 1069 
have unique kinematic distances and good fits. 

The luminosity distributions of these sources in terms of 
their RMS classification have been presented, and show that 
there are few MYSOs with L > 1O 5 L , though we detect UCH n 
regions up to ~1 x 1O 5 L , a difference which we find to be sta- 
tistically significant. We also find that using the SED fitter to 
obtain bolometric flux measurements is more consistent and has 
less scatter, than either using simple trapezium rule integration or 
greybody fits for mid-IR bright young massive stars. Finally, we 
obtained the flux ratios consistent with our SED fits, which al- 
low us to estimate the bolometric flux, of 280 sources for which 
far-IR fluxes could not be obtained. 

Having obtained the luminosity distributions of YSOs and 
H n regions in the RMS sample, we can go on to calculate the 
luminosity functions for these sources. This will be discussed in 
a forthcoming paper. 
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